package align_with_bl2seq;

require Exporter;


our @ISA = qw(Exporter);
our $VERSION = 1.00;
our @EXPORT = qw(	get_alignment_data
					expand_sequence_to_alignment
					expand_cDNA_with_alignment_array
					pos_in_aligned_seq
					expand_cDNA_to_alignment
					get_alignment_marks);

use strict;
use warnings;
use IO::File;
use POSIX qw(tmpnam);

# function declarations
sub parse_blast_results(\@);


my $tmp_file_name1;
my $tmp_file_name2;

sub delete_tmp_files
{
	if ($tmp_file_name1)
	{
		unlink $tmp_file_name1 or die "Couldn't unlink $tmp_file_name1 : $!";
	}
	if ($tmp_file_name2)
	{
		unlink $tmp_file_name2 or die "Couldn't unlink $tmp_file_name2 : $!";
	}
	$tmp_file_name1 = undef;
	$tmp_file_name2 = undef;

}
# delete files automatically when we exit or die
END
{
	delete_tmp_files();
}	


#_________________________________________________________________________________________

#	get_alignment_data

#		input:

#			$blast_prog			name and path of the blast programme
#			$seq1				sequence in fasta format
#			$seq2				sequence in fasta format

#		returns:

#			$score				bitscore
#			$alignment_length	length of the alignment
#			$e-value			e-value
#			$i_pos				reference to an array of position / length pairs
#			$j_pos				reference to an array of position / length pairs

#_________________________________________________________________________________________
sub get_alignment_data($\@\@)
{
	my ($blast_prog, $seq1, $seq2) = @_;

	# open temp files
	my ($tmp_fh1, $tmp_fh2);
	do {$tmp_file_name1 = tmpnam()}
		until $tmp_fh1 = IO::File->new($tmp_file_name1, O_RDWR | O_CREAT | O_EXCL );
	do {$tmp_file_name2 = tmpnam()}
		until $tmp_fh2 = IO::File->new($tmp_file_name2, O_RDWR | O_CREAT | O_EXCL );

	# write sequences to temp files and close
	print $tmp_fh1 @$seq1;
	print $tmp_fh2 @$seq2;

	close $tmp_fh1;
	close $tmp_fh2;

	# blast away
	
	my @results = split/\n/, `$blast_prog -i $tmp_file_name1 -j $tmp_file_name2 -p blastp -F f`;
	
	# can remove temp files now
	delete_tmp_files();

	# analyse blast results
	return (parse_blast_results(@results));

}

#_________________________________________________________________________________________

#	parse_blast_results

#		input:

#			@results
#				array of lines with the blast results

#		returns:

#			$score				bitscore
#			$alignment_length	length of the alignment
#			$i_pos				reference to an array of position / length pairs
#			$j_pos				reference to an array of position / length pairs

#_________________________________________________________________________________________
sub throw_bad_data($\@)
{
   	die "Error!\nUnexpected data on line ".$_[0]." of the blast output\n"
			.$_[1][$_[0]]."\n";
}
sub inc_line_with_check($\@)
{
	my $results = $_[1];
	DO_LOOP:
	{
		if (++$_[0] >= @$results)
		{
			die "Error!\nPremature end after line $[0]";
		}
		redo DO_LOOP if ($$results[$_[0]] eq "");
	}
}
sub parse_blast_results(\@)
{
	my ($results) = @_;
	my $score;
	my ($e_value, $len_a, $len_b);
	my ($alignment_length, $alignment_start);
	my (@i_pos, @j_pos);
	for (my $i = 0; $i < @$results; ++$i)
	{
		# get bit score		Score =  703 bits (1814), Expect = 0.0
		next unless ($results->[$i] =~ /^\s*Score\s*=\s*(\d+(?:\.\d+)?).*Expect\s*=\s*([0-9\-\.Ee]+)/);
		$score = $1;
		$e_value = $2;
		inc_line_with_check($i, @$results);
	
		# get the length of the alignment i.e. "Identities = 527/625 (84%)"
			throw_bad_data($i, @$results) unless
		$results->[$i] =~ /^\s*Identities\s*=\s*\d+\/(\d+)/;
		$alignment_length = $1;
		inc_line_with_check($i, @$results);
	
		# go through query/align/subject line triplets of the top alignment
		my $do_initial_align = 1;
		while ($results->[$i] =~ /^\s*Query:/)
		{
			# query
				throw_bad_data($i, @$results) unless
			$results->[$i] =~ /^\s*					# whitespace
								Query:				# name of line
								\s*                 # whitespace
								(\d+)               # start position of query
								\s*                 # whitespace
								([A-Z\-+]+)/x;		# text string
	
			my $i_aa_line	= $2;
			my $i_seq_start = $1;
			inc_line_with_check($i, @$results);
	
			#skip past alignment line
			inc_line_with_check($i, @$results);
	
	
			# subject
			throw_bad_data($i, @$results) unless
				$results->[$i] =~ /^\s*					# whitespace
									Sbjct:				# name of line
									\s*                 # whitespace
									(\d+)               # start position of query
									\s*                 # whitespace
									([A-Z\-+]+)/x;		# text string
			my $j_aa_line = $2;
			my $j_seq_start = $1;
			inc_line_with_check($i, @$results);
	
			# note the offset to the first aligned position
			if ($do_initial_align)
			{
	
				$alignment_start = $i_seq_start;
				if ($i_seq_start > $j_seq_start)
				{
					push(@j_pos, [0, $i_seq_start-$j_seq_start]);
				}
				elsif ($i_seq_start < $j_seq_start)
				{
					push(@i_pos, [0, $j_seq_start-$i_seq_start]);
					$alignment_start = $j_seq_start;
				}
				$do_initial_align = 0;
			}
	

			# save the number and position of the aligning characters for each sequence
			while ($i_aa_line =~ /(\-+)/g)
			{
				# save position and length of aligning markers
				my $len = length($1);
				my $pos = pos($i_aa_line) - $len + $i_seq_start - 1;
				if (@i_pos && $i_pos[-1][0] == $pos)
				{
					$i_pos[-1][1] += $len;
				}
				else
				{
					push(@i_pos, [$pos, $len]);
				}
				# adjust starting position to account for the markers already counted
				$i_seq_start -= $len;
			}
	
			while ($j_aa_line =~ /(\-+)/g)
			{
				# save position and length of aligning markers
				my $len = length($1);
				my $pos = pos($j_aa_line) - $len + $j_seq_start - 1;
				if (@j_pos && $j_pos[-1][0] == $pos)
				{
					$j_pos[-1][1] += $len;
				}
				else
				{
					push(@j_pos, [$pos, $len]);
				}
				# adjust starting position to account for the markers already counted
				$j_seq_start -= $len;
			}
		}
		
		if (defined $score)
		{
			#find length of two sequences
			inc_line_with_check($i, @$results)
				while ($results->[$i] !~ /^length of query:\s*(\d+)/);
			$len_a = $1;
			inc_line_with_check($i, @$results)
				while ($results->[$i] !~ /^length of database:\s*(\d+)/);
			$len_b = $1;
		}
		my $final_len_a = $len_a;
		$final_len_a += $_->[1] for (@i_pos);
		my $final_len_b = $len_b;
		$final_len_b += $_->[1] for (@j_pos);

		# add end alignment marks
		if ($final_len_a > $final_len_b)
		{
			push(@j_pos, [$len_b, $final_len_a - $final_len_b]);
		}
		elsif ($final_len_b > $final_len_a)
		{
			push(@i_pos, [$len_a, $final_len_b - $final_len_a]);
		}
=pod
		#debug		
		print STDERR "$_->[0],$_->[1];" for (@i_pos);
		print STDERR "\n";
		print STDERR "$_->[0],$_->[1];" for (@j_pos);
		print STDERR "\n";
=cut
		return ($score, $alignment_start, $alignment_length, $e_value, $len_a, $len_b, \@i_pos, \@j_pos,
					$results);
	}
	return;

}




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	print the modified sequences

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#_________________________________________________________________________________________

#	expand_cDNA_to_alignment

#		input:

#			$seq
#			$aligning character
#			@list of position/length pairs in string format 0:63|119:6|263:1

#		returns:

#			$expanded sequences
#_________________________________________________________________________________________
sub expand_cDNA_with_alignment_array($$$)
{
	my ($seq, $char, $align) = @_;
	$char |= '-';
	$char = substr($char, 0,  1);
	$char = $char x 3;
	
	foreach (reverse @$align)
	{
		substr($seq, $_->[0] * 3, 0) = ($char x $_->[1]);
	}
	return $seq;

}

sub expand_cDNA_to_alignment($$$)
{
	my ($seq, $char, $align_str) = @_;
	my @align = map {[split '\:', $_]; } split '\|', $align_str;
	return expand_cDNA_with_alignment_array($seq, $char, @align)
}


#_________________________________________________________________________________________

#	expand_sequence_to_alignment

#		input:

#			$seq
#			$aligning character
#			@list of position/length pairs

#		returns:

#			$expanded sequences
#_________________________________________________________________________________________
sub expand_sequence_to_alignment($$\@)
{
	my ($seq, $char, $align) = @_;
	$char |= '-';
	$char = substr($char, 0,  1);
	
=pod
		#debug		
		print STDERR "length = ", length $seq, "\n";
		print STDERR "$_->[0],$_->[1];" for (@$align);
		print STDERR "\n";
=cut
	foreach (reverse @$align)
	{
		substr($seq, $_->[0], 0) = ($char x $_->[1]);
	}
	return $seq;

}


#_________________________________________________________________________________________

#	pos_in_aligned_seq

#		input:

#			$position
#			@list of alignment position/length pairs

#		returns:

#			adjusted position given the aligning marks
#_________________________________________________________________________________________
sub pos_in_aligned_seq($\@)
{
	my ($pos, $align) = @_;
	my $new_pos = $pos;
	foreach (@$align)
	{
		$new_pos += $_->[1] if ($pos >= $_->[0]);
	}
	return $new_pos;

}


#_________________________________________________________________________________________

#	pos_from_aligned_seq

#		input:

#			$position
#			@list of alignment position/length pairs

#		returns:

#			adjusted position given the aligning marks
#_________________________________________________________________________________________
sub pos_from_aligned_seq($\@)
{
	my ($pos, $align) = @_;
	foreach (@$align)
	{
		# before next gap
		last if ($_->[0] >= $pos);
		
		$pos -= $_->[1];

		# in gap region
		if ($_->[0] >= $pos)
		{
			$pos = $_->[0];
			last;
		}
	}
	return $pos;
}

#_________________________________________________________________________________________

#	pos_via_alignment


#		map from one alignment to another


#		input:

#			$position
#			@list of alignment position/length pairs for my sequence
#			@list of alignment position/length pairs for other sequence

#		returns:

#			adjusted position given the aligning marks
#_________________________________________________________________________________________
sub pos_via_alignment($\@\@)
{
	my ($pos,  $align1,  $align2) = @_
	return pos_from_aligned_seq(pos_in_aligned_seq($pos, @$align1),  @$align2));

}

#_________________________________________________________________________________________

#	get_alignment_marks

#		input:

#			$ aligned protein sequence with no spaces etc.

#		returns:

#			positions of aligning marks
#_________________________________________________________________________________________
sub get_alignment_marks($)
{
	my ($seq) = @_;
	my @align_pos;
	my $cummulative = 0;
	# save the number and position of the aligning characters for each sequence
	while ($seq =~ /(\-+)/g)
	{
		# save position and length of aligning markers
		my $len = length($1);
		push(@align_pos, [pos($seq) - $len - $cummulative, $len]);
		$cummulative += $len;
	}
	return @align_pos
}

1;

